---
title: "Analytical Paleobiology Tutorials 2018"
subtitle: "Diversity analyses - richness, diversity, evenness"
author: Jessica Blois
date: "August 1, 2018 (last updated: October 29, 2018)"
output:
  pdf_document: default
  html_notebook: default
---
# Part 1: Introduction to estimating biological diversity

### What is the goal?  

Often, we are interested in describing the *richness* and structure or *diversity* of an assemblage. But more importantly, we are often more interested in comparing richness and diversity across space and/or time, by merging data from multiple sites or from different sampling units within a site.  

By *richness*, we mean the number of unique species or taxonomic/phylogenetic/functional units in the assemblage.  By *diversity*, we usually mean some aspect of the structure of the assemblage - a combination of both the number of unique elements as well as the relative abundance of the different elements. 

Overall, biological diversity metrics can be characterized as:

1) Richness metrics : number of species
2) Evenness metrics : capture the distribution of abundance among species
3) Diversity metrics : combines both richness and evenness 
4) Dominance metrics : somewhat the inverse of evenness, focused on the most common species (not covered)

Today, we will be exploring different aspects of estimating richness and diversity. A good resource for this topic, and one on which I have relied heavily in this lecture (in addition to package vignettes), is the book "Biological Diversity" by Magurran and McGill <https://global.oup.com/academic/product/biological-diversity-9780199580675?cc=us&lang=en&>.  This book may be available through your institutional library.

# Part 2: Getting started with R

We will need to load all relevant package 
```{r}
# uncomment the following line if you haven't installed these packages yet
# install.packages(c('vegan', 'neotoma', 'dplyr', 'sads', 'iNEXT')) 

library('vegan')
library('neotoma')
library('dplyr')
library('sads')
library('iNEXT')
```

# Part 3: Diversity analyses for a single site through time

### Download data for a single site

Let's start with a simple search for a single site in the Neotoma Paleoecology Database <https://www.neotomadb.org/>. We'll download mammal data from one site (late Quaternary cave deposit) in northern California, Samwell Cave Popcorn Dome, along with the corresponding ages of the sampled units.
```{r}
samwell_site <- get_site(sitename="Samwell Cave") #there's one site in the database that matches this
samwell_datasets <- get_dataset(samwell_site[[1]]) # this search lists all datasets associated with Samwell Cave. There are 7 different datasets for the site: 4 containing vertebrate faunal information and 3 containing geochronologic information. 
browse(samwell_datasets) # Note that this opens a new window in your browser

# Download info for just the last dataset, 14262. 
scpd_dataset <- get_dataset(samwell_site[[1]])$'14262'
scpd_download <- get_download(scpd_dataset)

# Extract the count data and the sample ages
scpdFauna <- counts(scpd_download)[[1]]
scpdAges <- ages(scpd_download)[[1]]
   
# save as R object for later upload in case no internet or API not working
# save(scpdFauna, scpdAges, file="data/scpdObjects.RData")

# Note: Dataset 14261 contains the underlying radiocarbon dates for the site, in case you want to examine those. You'll need to uncomment the next two lines
# scpd_geochron <- get_geochron(14261)
# scpd_geochron[[1]]$geochron
```

### Clean up the assemblage/site data

OK, now that we have the basic dataset downloaded, we will need to clean up the data. 

**Question:**

What kinds of data cleaning might we want to do here?

It's often helpful to examine the structure of the data to answer this question.
```{r}
# examine the final downloaded file
str(scpdFauna) # another possible way to examine the data: glimpse(scpdFauna)
head(scpdFauna) 
```

Start cleaning up the data files.
```{r}
# What are the rownames showing?  
# Let's replace the sample ID in the fauna object (rownames) with the ages
# but only do this if the sample IDs match across objects
ages <- scpdAges$age
if (all(rownames(scpdFauna) == scpdAges$sample.id)){
  rownames(scpdFauna) <- ages
}

# set the names of the small mammal genera of interest
# note: this is a clunky workaround - ideally we'd be able to harmonize taxonomy using the database/API tools, but they aren't as good in Neotoma as in PBDB because of the way taxonomy is handled
sm_genus <- c("Scapanus", "Microtus", "Neotoma", "Peromyscus", "Sorex", "Spermophilus", "Tamias", "Tamiasciurus", "Thomomys", "Arborimus", "Reithrodontomys", "Sciurus", "Glaucomys", "Myodes", "Aplodontia")

# create a new dataframe with the merged data
scpdSM <- matrix(nrow=nrow(scpdFauna), ncol=length(sm_genus))

for (i in 1:length(sm_genus)){
  cols <- grep(sm_genus[i], colnames(scpdFauna))

  if (length(cols)==1){
    scpdSM[,i] <- scpdFauna[,cols]
  }
  if (length(cols) > 1){
    scpdSM[,i] <- rowSums(scpdFauna[,cols])
  }
  
}

# clean up dataframe
colnames(scpdSM) <- sm_genus
rownames(scpdSM) <- rownames(scpdFauna)
save(scpdSM, file="data/scpdFauna_cleaned.RData")
```

### Richness analyses

**Rarefaction-based estimates**

First, let's calculate raw richness, for each sample through time
```{r}
# we can practice building our own functions to count up the species
countFcn <- function(df){
  count <- length(which(df>0))
  return(count)
} 
raw.S <- apply(scpdSM, 1, countFcn)
# or you can do the same calculation using a built-in function (specnumber) in vegan
specnumber(scpdSM)
```

However, we know these are fossil data and (even if they weren't!) sampling effort may have varied through time and across samples. ie, the total sample size (N) of each sample (stratum) varies substantially, as can be seen by examining the range of total N.
```{r}
range(rowSums(scpdSM))
minN <- min(rowSums(scpdSM))
```

So one approach is rarefy each assemblage down to the minimum N when calculating richness. This effectively answers: What is the expected richness of an assemblage if we sampled the same number of individuals as were observed in some other assemblage?  You did this manually in Seth's section, and then he also introduced you to the function *rarefy* in the *vegan* package. *rarefy* "gives the expected species richness in random subsamples of size sample from the community" (from help doc) and performs this for all samples/assemblages in the community matrix.
```{r}
richness <- rarefy(scpdSM, sample=minN, se=T) # because I specified *se=T*, this function outputs both S and se
S <- richness[1,]
S.se <- richness[2,]

# plot the result through time using base R plotting function
plot(S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, floor(max(S+S.se))+1), xlab="Age (cal BP)", ylab=expression('Richness' %+-% 'se'))
#arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)
suppressWarnings(arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)) 

# Note: use suppressWarnings with caution.  I added the suppressWarnings tag so that the only thing output is the plot. But I only did this after I confirmed that the warning message was ok. Switch the arrows lines that are commented/uncommented and re-run the plot with that line of code for arrows to view the warning and decide for yourself whether it should be suppressed or not.
```

Let's compare rarefied to raw richness, for each sample through time
```{r}
plot(raw.S ~ ages, type="l", xlim=c(21000, 0), ylim=c(0, max(raw.S)+1), lty=2, col="gray", xlab="Age (cal BP)", ylab='Richness')
lines(S ~ ages, type="l")
arrows(ages, S+S.se, ages, S-S.se, code=3, length=0.1, angle=90)
```

We can also compare the sampling curves using function *rarecurve*. This grabs random samples of different numbers of individuals in a step-wise function, then calculates richness for each sample.
```{r}
# set colors, from cool for the older (LGM) times and red for the warmer (Holocene) times
ageColors <- rainbow(length(ages), start=0, end=4/6)

rarecurve(scpdSM, col=ageColors)
```

**Extrapolation-based estimates of richness**

Rarefaction calculates the expected richness if the same number of individuals were sampled across different samples. However, there are other ways to estimate richness and compare among samples. Among the more common alternatives to rarefaction are the Chao and ACE (Abundance-based Coverage Estimator) estimators of richness. Let's first calculate them, and then discuss.

```{r}
?estimateR # read more about this function, and how ACE and Chao1 are calculated
alt_richness <- estimateR(scpdSM)

# Bind all richness estimates together and compare:
richness_all <- rbind(alt_richness, richness)

# plot them against one another
y.max <- floor(max(richness_all))+1 # set y max for plotting 

mainEstimates <- c('S.obs', 'S.chao1', 'S.ACE', 'S') 
plotRows <- match(mainEstimates, rownames(richness_all))
plotCols <- c("black", "blue", "green", "red")
  
i=1
plot(richness_all[plotRows[i],] ~ ages, type="l", ylim=c(0, y.max), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[i])

for (i in 2:length(plotRows)){
    lines(richness_all[plotRows[i],] ~ ages, col=plotCols[i])
}
legend("bottomleft", legend=mainEstimates, lty=1, col=plotCols, bty="n", cex=0.75)
```

**Discussion**

You can see that these two approaches (rarefaction vs extrapolation) produce fundamentally different estimates of the trend of species richness through time.  Individual-based rarefaction (S) shows a trend towards lower richness from the Pleistocene to the Holocene, whereas ACE and Chao1 estimates show a trend towards higher richness from the Pleistocene to the Holocene.  

So what's the difference?  Individual-based rarefaction (S) estimates expected richness given a subsample of n individuals from the assemblage. It's a simple form of standardization so that different samples are comparable. In this case, n is sampled at the minimum N observed among the samples (N=130).  And since n is generally lower than the total N sampled, it's really an underestimate of on-the-ground richness.

On the other hand, Chao1 and ACE estimate the asymptotic richness estimates. ie, they both try to extrapolate beyond the sampled individuals to figure out the expected richness if the assemblage was sampled completely. They focus in particular on information in the detected rare species (singletons, etc) to estimate undetected species. 

Chao1 and ACE have been shown to perform better than other estimators in various studies. They are also similar to (or fundamentally the same as) SQS. A recent study (Close et al. 2018) examined differences among richness estimators and found "sampling in the tetrapod fossil record is generally not yet good enough to use extrapolators unless they are applied within a rarefy‐and‐extrapolate protocol of the kind used in our simulation experiments" and "In particular, rarefying extrapolators to equal coverage is the best way to remove confounding effects of among‐sample variation in evenness, a problem that affects all richness estimators when sampling is comparatively limited". The package *iNEXT* implements this rarefying extrapolators for Chao1, so let's do a quick tour through that package.

```{r}
# this chunk uses the iNext package
vignette("Introduction", package="iNEXT")

iNextRichness <- iNEXT(t(scpdSM)) 
# Note 1: we had to transpose the community matrix
# Note 2: this calculates not just richness, but also Shannon and Simpson diversities.

# You can also rarefy within iNEXT using estimateD
D_minN <- estimateD(t(scpdSM), datatype="abundance", base="size", level=NULL) 
D_minN
# Note 1: SC = sample coverage, qD is the estimate with lower and upper confidence levels
# Note 2: Richness estimates here are the same as using *rarefy* in vegan. e.g., compare the following:
D_minN[which(D_minN$order == 0),]$qD
richness[1,]
```

Plot the Chao estimators
```{r}
ggiNEXT(iNextRichness)
ggiNEXT(iNextRichness, type=1, facet.var="site")
ggiNEXT(iNextRichness, type=2)
```

### Other diversity measures

OK, up to now we've just been counting the number of species in different ways.  The other aspect of assemblages is their abundance structure - how abundance varies among the different species in the deposit.

**Diversity**

There are several common diversity measures, which combine both richness and evenness.  I think it's helpful to view the calculations for these by looking at one of the vegan vignettes
```{r}
vignette("diversity-vegan") # note: this opens a pdf

H <- diversity(scpdSM, index="shannon")
Simp <- diversity(scpdSM, index="simpson")
```

Practice calculating Shannon on your own
```{r}
# Shannon
pi.SCPD <- scpdSM/rowSums(scpdSM)
pi.SCPDxlogpi <- pi.SCPD * log(pi.SCPD)
sums_shan <- rowSums(pi.SCPDxlogpi, na.rm=T)
H_calc <- -sums_shan
all(H==H_calc)

# Simpson
pi.SCPD2 <- pi.SCPD^2
sums_simp <- rowSums(pi.SCPD2, na.rm=T)
Simp_calc <- 1-sums_simp
all(Simp==Simp_calc)
```

**Evenness**

You were introduced to evenness is Seth's lecture. There aren't good built-in functions for calculating evenness (that I know of), but Pielou's (or Shannon) evenness can be derived from Shannon diversity (H). (H/ln(S)).  The idea underlying both of these measures is that if diversity incorporates both richness and evenness, then factoring out richness from the diversity measures should leave us with evenness.
```{r}
Shan_even <- H/log(specnumber(scpdSM))
Simp_even <- Simp/log(specnumber(scpdSM))
```

Let's plot these through time
```{r}
plotCols <- c("black", "blue",  "green", "red")
plot(H~ages, type="l", ylim=c(0, 2), xlim=c(21000, 0), xlab="Age (cal BP)", ylab="Richness", col=plotCols[1])
lines(Simp~ages, col=plotCols[2])
lines(Shan_even~ages, col=plotCols[3], lty=2)
lines(Simp_even~ages, col=plotCols[4], lty=2)
legend("bottomleft", legend=c("Shannon", "Simpson", "Evenness_Shannon", "Evenness_Simpson"), lty=c(1,1,2,2), col=plotCols, bty="n", cex=0.75)
```

Note: you could also calculate Shannon and Simpson through iNext (as above), and probably using several other packages.


# Part 4: Extending these analyses to multiple sites
Let's ask: How do richness, evenness, and diversity vary across latitude and longitude? We will use a dataset of Quaternary fossil pollen data from eastern North America for 6000 years before present.

Note, this is not count data - the data are provided as the relative abundance of each taxon at the site - so the rarefaction and extrapolation estimators don't work. They only accept count data / integers.

**Data loading and prep**
```{r}

# load the data and examine the structure
pollen.6000 <- read.csv(file="data/pollen/6000.interp.data_clean_sites.csv")
str(pollen.6000)

# clean up the site/rownames
pollen.6000.sites <- pollen.6000$sites
rownames(pollen.6000) <- pollen.6000$sites
pollen.6000 <- pollen.6000[,-1]
```

**Richness, diversity, and evenness calculations**
```{r}
# calculate richness and diversity
S.obs <- specnumber(pollen.6000)
H <- diversity(pollen.6000, "shannon")
Simp <- diversity(pollen.6000, 'simpson')
H_even <- H/log(S.obs)
```

OK, now examine how diversity changes across lat and long
```{r}
# get site lat/longs
meta <- read.delim(file="data/pollen/All.site.data-withagemodel-finalv2.txt", sep="\t")
site.longlats <- meta[match(pollen.6000.sites, meta$Handle), c('Longitude', 'Latitude')]


# simple relationships with latitude and longitude
par(mfrow=c(2,2), mar=c(4,4,1,1)+0.1)
plot(H~site.longlats$Latitude, pch=16, xlab="Latitude")
plot(S.obs~site.longlats$Latitude, pch=16, xlab="Latitude")
plot(Simp~site.longlats$Latitude, pch=16, xlab="Latitude")
plot(H_even~site.longlats$Latitude, pch=16, xlab="Latitude")

par(mfrow=c(2,2), mar=c(4,4,1,1)+0.1)
plot(H~site.longlats$Longitude, pch=16, xlab="Longitude")
plot(S.obs~site.longlats$Longitude, pch=16, xlab="Longitude")
plot(Simp~site.longlats$Longitude, pch=16, xlab="Longitude")
plot(H_even~site.longlats$Longitude, pch=16, xlab="Longitude")

```

### Practice
Note: In the 'data' folder, you have similar pollen files for multiple time periods.

**Questions:** 

1. What else could you do with these data? Explore additional diversity analyses, either using these data for 6000 years BP, or by looking at the datasets for other times.

2. Use the knowledge you've gained from today or previous modules in the course to extract information from other sites in the PaleobioDB or Neotoma and calculate diversity.

# Part 5: SADs and RADs
One other common way to visualize the structure of communities is through species abundance distributions (SADs), and related rank abundance distributions (RADs). This is a way to visualize the differences in commonness and rarity among species, and can potentially be compared among communities. SADs/RADs have a long history in ecology, and the general pattern is that to be rare is common and to be common is rare.  That is, very few species are highly abundant, and many more species have low abundances.

```{r}
# this chunk uses the sads package
vignette('sads_intro') # note: opens a pdf
```

Let's revisit the small mammal data from Samwell Cave and answer the question: has the structure of the community changed through time? We'll focus mostly on RADs.
```{r}

# examine the two endpoints
scpd1.oc <- octav(scpdSM[1,])
scpd14.oc <- octav(scpdSM[14,])
scpd1.oc
scpd14.oc

par(mfrow=c(4,4), mar=c(3,3,2,0)+0.1)
for (i in 1:nrow(scpdSM)){
  plot(rad(scpdSM[i,]), prop=F, xlab="", ylab="", xlim=c(1, 13), type="b", pch=16)
  mtext(side=2, line=2, "Number of Indiv", cex=0.5)
  mtext(side=3, line=0.5, paste0("Time=", ages[i]), cex=0.75)
  mtext(side=1, line=2, "Rank", cex=0.5)
}
```

**Question:**

What do the figures above show us?

**Alternate approach**

McGill (2011, Ch. 9) identifies a few problems with simplistic visual comparisons of RADs that are visible in the figure above.  

1) Ranks can be difficult to deal with mathematically
2) Because richness varies among communities, it can be difficult to compare RADs across assemblages. It's not such a huge deal with the Samwell Cave data, but would be a bigger issue with more disparity in richness.

He suggests instead using ECDFs (empirical cumulative density functions) instead.
```{r}
par(mfrow=c(4,4), mar=c(3,3,2,0)+0.1)
for (i in 1:nrow(scpdSM)){
  plot(ecdf(scpdSM[i,]), xlab="", ylab="", xlim=c(1, 13), main="", verticals=T, do.points=F, cex.axis=0.75)
  mtext(side=2, line=2, "% of Species", cex=0.5)
  mtext(side=3, line=0.5, paste0("Time=", ages[i]), cex=0.75)
   mtext(side=1, line=2, "% abundance", cex=0.5)
}
```

**Formal model fitting**

It's possible using the *sads* package (and also *vegan*, though the approaches differ) to fit various models to the SADs/RADs. There are quite a few different models to choose from, and different authors have postulated that different distributions are indicative of different types of communities. I've randomly chosen three: geometric (low diversity), Fisher's log-series (medium diversity), and lognormal (high diversity).
```{r}
# fit the models
scpd1.geom <- fitsad(scpdSM[1,scpdSM[1,]>0], 'geom')
scpd1.ls <- fitsad(scpdSM[1,scpdSM[1,]>0], 'ls')
scpd1.lnorm <- fitsad(scpdSM[1,scpdSM[1,]>0], 'lnorm')

# compare the fit of each model
AICtab(scpd1.geom, scpd1.ls, scpd1.lnorm, base=TRUE)

# transform to appropriate class for plotting
scpd1.geom.oc <- octavpred(scpd1.geom)
scpd1.ls.oc <- octavpred(scpd1.ls)
scpd1.lnorm.oc <- octavpred(scpd1.lnorm)

# plot the three models against one another and the data
plot(scpd1.oc)
lines(scpd1.geom.oc, col="red")
lines(scpd1.ls.oc, col="blue")
lines(scpd1.lnorm.oc, col="green")
legend("topright", c("geometric", "Logseries", "Lognormal"), lty=1, col=c("red","blue", "green"))
```
